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The hybrid kinetic model supports comprehensive simulation of the interaction between different 
spatial and energetic elements of the Europa moon-magnetosphere system with respect to a variable 
upstream magnetic field and flux or density distributions of plasma and energetic ions, electrons, and 
neutral atoms. This capability is critical for improving the interpretation of the existing Europa flyby 
measurements from the Galileo Orbiter mission, and for planning flyby and orbital measurements 
(including the surface and atmospheric compositions) for future missions. The simulations are based on 
recent models of the atmosphere of Europa (Cassidy et al., 2007; Shematovich et al., 2005). In contrast 
to previous approaches with MHD simulations, the hybrid model allows us to fully take into account 
the finite gyroradius effect and electron pressure, and to correctly estimate the ion velocity distribution 
and the fluxes along the magnetic field (assuming an initial Maxwellian velocity distribution for 
upstream background ions). Photoionization, electron-impact ionization, charge exchange and colli- 
sions between the ions and neutrals are also included in our model. We consider the models with 0 + + 
and S ++ background plasma, and various betas for background ions and electrons, and pickup 
electrons. The majority of 0 2 atmosphere is thermal with an extended non-thermal population 
(Cassidy et al., 2007). In this paper, we discuss two tasks: (1) the plasma wake structure dependence 
on the parameters of the upstream plasma and Europa’s atmosphere (model I, cases (a) and (b) with a 
homogeneous Jovian magnetosphere field, an inductive magnetic dipole and high oceanic shell 
conductivity); and (2) estimation of the possible effect of an induced magnetic field arising from 
oceanic shell conductivity. This effect was estimated based on the difference between the observed and 
modeled magnetic fields (model II, case (c) with an inhomogeneous Jovian magnetosphere field, an 
inductive magnetic dipole and low oceanic shell conductivity). 

© 2013 Elsevier Ltd. All rights reserved. 


1. Introduction 

The interaction of the Jovian plasma torus with Europa and 
other moons is a fundamental problem in magnetospheric physics 
(see e.g., Goertz, 1980; Southwood et al., 1980; Southwood and 
Dunlop, 1984; Wolf-Gladrow et al., 1987; Ip, 1990; Schreier et al., 
1993; Lellouch, 1996). The plasma environment near Europa was 
studied by flyby observations during the Galileo prime mission 
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and the extended Galileo Europa mission (Kivelson et al., 1997; 
Khurana et al., 1998; Kivelson et al., 1999; Paterson et al., 1999). 

Europa, one of the icy moons of Jupiter, was encountered by 
the Galileo satellite three times during its primary mission, seven 
times during its Galileo Europa Mission (GEM), and once during 
Galileo Millennium Mission (GMM). Europa is located at a radial 
distance of 9.4R; (Jovian radii, 71,492 km) from Jupiter and has a 
radius of 1560 km (1 R E ). 

The interaction of Europa with the magnetized plasma of the 
Jovian plasma sheet gives rise to a so-called Alfven wing, which 
has been extensively studied in the case of Io (e.g., Neubauer, 
1980; Southwood et al., 1980; Herbert, 1985; Lipatov and Combi, 
2006). Neubauer (1998, 1999) has shown theoretically how an 
Alfven wing is modified by an induced magnetic field, such as that 
found at Europa (Kivelson et al., 2000). Observations by Kivelson 
et al. (1992) show the generation of ultra-low frequency electro- 
magnetic waves in Europa’s wake. These waves have frequencies 
near and below the gyrofrequencies of the ion species in the 
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plasma torus (e.g M ionized sulfur, oxygen, and protons). Ion 
cyclotron waves grow when ion distribution functions are suffi- 
ciently anisotropic, as occurs when ion pickup creates a ring 
distribution of ions (in velocity space). The analysis of these 
waves has been done by Huddleston et al. (1997) (Io), Volwerk 
et al. (2001) and Kivelson et al. (2009) (Europa). They found 
intensive wave power at low frequencies (near and below the 
cyclotron frequencies of heavy ions) in Europa’s wake during the 
Ell and El 5 flybys. However, our current 3D hybrid modeling 
cannot yet produce these waves due to insufficient spatial grid 
resolution. 

The most general and accurate theoretical approach to this 
problem would require the solution of a nonlinear coupled set of 
integro-MHD/kinetic-Boltzmann equations which describe the 
dynamics of Jupiter’s corotating magnetospheric plasma, pickup 
ions, and ionosphere, together with the neutrals from Europa’s 
atmosphere. To first order, the plasma and neutral atoms and 
molecules are coupled by charge exchange and ionization. The 
characteristic scale of the ionized components is usually deter- 
mined by the typical ion gyroradius, which for Europa is much 
less than characteristic global magnetospheric scales of interest, 
but which may be comparable to the thickness of the plasma 
structures near Europa. Kinetic approaches, such as Direct Simu- 
lation Monte Carlo, have been applied to the understanding of 
global aspects of the neutral atmosphere (Marconi et al., 1996; 
Austin and Goldstein, 2000). Plasma kinetic modeling is, however, 
much more complicated, and even at the current stage of 
computational technology requires some approximations and 
compromises to make some initial progress. Several approaches 
have been formulated for including the neutral component and 
pickup ions self-consistently in models that describe the interaction 
of the plasma torus with Europa. 

There have been recent efforts to improve and extend the 
pre-Galileo models for Europa, Io and Ganymede, in terms of the 
MHD (Kabin et al., 1999; Combi et al., 1998; Linker et al., 1998; 
Kabin et al., 2001; Jia et al., 2008), the electrodynamic (Saur et al., 
1998, 1999; Schilling et al., 2008), and hybrid kinetic (Lipatov and 
Combi, 2006; Lipatov et al., 2010) approaches. These approaches 
are distinguished by the physical assumptions that they include. 
MHD and hybrid kinetic models cannot, at least yet, include the 
charge separation effects which are likely to be important very 
close to the moon where the neutral densities are large and the 
electric potential can introduce non-symmetric flow around the 
body. MHD models for Io either include constant artificial con- 
ductivity (Linker et al., 1998) or assume perfect conductivity 
(Combi et al., 1998). Comparisons of the sets of published results 
do not indicate that this choice has any important consequences. 
The MHD model of Europa developed by Kabin et al. (1999) 
includes an exospheric mass loading, ion-neutral charge exchange, 
and recombination. Further development of this model by Liu 
et al. (2000) already includes a possible intrinsic dipole magnetic 
field of Europa. Schilling et al. (2007, 2008) found that for the 
conductivity of Europa’s ocean values of 500 mS/m or large 
combined with ocean thickness of 100 km and smaller to be most 
suitable to explain the magnetic flyby data. They also found that 
the influence of the fields induced by the time variable plasma 
interaction is small compared to the induction caused by the 
time-varying background field. 

Hybrid kinetic models can include the finite ion gyroradius 
effects, non-Maxwellian velocity distribution for ions, and correct 
flux of pickup ions along the magnetic field. Hybrid modeling of Io 
has demonstrated several features. The kinetic behavior of ion 
dynamics reproduces the inverse structure of the magnetic field 
(due to drift current) which cannot be explained by standard 
MHD or electrodynamic modeling which do not account for 
anisotropic ion pressure. The diamagnetic effect of non-isotropic 


gyrating pickup ions broadens the B-field perturbation and 
produces increased temperatures in the flanks of the wake, as 
observed by the Galileo spacecraft, but had not been explained by 
previous models. The temperatures of the electrons which are 
created and cooled by collisions with neutrals in the exosphere 
and inside the ionosphere may strongly affect the pickup ion 
dynamics along the magnetic field and consequently the pickup 
distribution across the wake. The physical chemistry in Io’s 
corona was considered in the paper by Dols et al. (2008). They 
couple a model of the plasma flow around Io plus a multi-species 
chemistry model and compare the model results to the Galileo 
observation in Io’s wake. 

Galileo flyby measurements E4, E6 (plasma only), Ell, El 2, 
El 4, El 5, El 9, and E26 demonstrate several features in the plasma 
environment: Alfven wing formation and an induced magneto- 
sphere, possible existence of the dipole-type induced magnetic 
field, and variation of the magnetic field in the plasma wake due 
to diamagnetic currents. The measurements also demonstrate 
mass loading of the plasma torus plasma by pickup ions and the 
interaction of the ions with the surface of Europa. For an 
interpretation of these data, we need to use a kinetic model 
because of effects of the finite ion gyroradius. 

Hybrid models have been shown to be very useful in studying 
the complex plasma wave processes of space, astrophysical, and 
laboratory plasmas. These models provide a kinetic description of 
plasmas in local regions, together with the possibility of perform- 
ing global modeling of the whole plasma system. Revolutionary 
advances in computational speed and memory are making hybrid 
modeling of various space plasma problems a much more effec- 
tive general tool. 

In this paper, we apply a time-dependent Boltzmann equa- 
tion (a “particle in cell” approach) together with a hybrid kinetic 
plasma (ion kinetic) model in three spatial dimensions (see, e.g. 
Lipatov and Combi, 2006; Lipatov et al., 2010), using a prescribed 
but adjustable neutral atmosphere model for Europa. A Boltz- 
mann simulation is applied to model charge exchange between 
incoming and pickup ions and the immobile atmospheric neu- 
trals. In this paper, we discuss the results of the hybrid kinetic 
modeling of Europa’s environment — namely the global plasma 
structures (formation of the magnetic barrier, Alfven wing, pickup 
ion tail, etc.). The results of these kinetic modeling are compared 
with the Galileo E4 flyby observational data. Currently, we are 
working on the hybrid model of the El 2 flyby. The remarkable 
aspect of this flyby is a strong variation in the upstream plasma 
density profile approximately from 400 cm -3 to 80 cm -3 . The 
results of this modeling will be discussed in future publications. 

The paper is organized as follows: in Section 2 we present 
the computational model and a formulation of the problem. 
In Section 3 we present the results of the modeling of the plasma 
environment near Europa and the comparison with observational 
data. Finally, in Section 4 we summarize our results and discuss 
the future development of our computational model. 


2. Formulation of the problem and mathematical model 

To study the interaction of the plasma torus with the ionized 
and neutral components of Europa’s environment, we use a 
quasineutral hybrid model for ions and electrons. The model 
includes ionization (which in the Europa environment is domi- 
nated by electron impact ionization, not photoionization) and 
charge exchange. The atmosphere is considered to be an immo- 
bile component in this paper. 

In our hybrid modeling, the dynamics of upstream ions 
and implanted ions are described in a kinetic approach, while 
the dynamics of the electrons are described in a hydrodynamical 
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approximation. The details of this plasma-neutral approach were 
developed early for the study of the Io-Jovian plasma interaction 
(Lipatov and Combi, 2006). 

The single ion particle motion is described by the equations 
(see, e.g. Eqs. (1) and (14) from Mankofsky et al., 1987) 


dt 


= v s>/ , 


dv S ,/ 

dt 


= W t (E+ 


v s>/ X B' 


m e vu 

' Mi 


(v,i-U,)-^J-v io v s ,, 
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( 1 ) 


Here we assume that the charge state is Z,=l. U z , and J denote the 
charge-averaged velocity of all (incoming and pickup) ions and 
the total current, Eq. (5). The subscript s denotes the ion popula- 
tion (s=l,2 for incoming ions and s = 3,4 for pickup ions) and the 
index l is the particle index. v ie and v io are collision frequencies 
between ions and electrons, and ions and neutrals that may 
include Coulomb collisions and collisions due to particle-wave 
interaction. 

For a plasma, the thermal velocity, v' a (a = i,e), is assumed 
greater than the drift velocity, so we take 

v a ,o = n 0 a°’ a v' a , (2) 

where the cross section er°’ a is typically about 5 x 10“ 15 cm 2 (see, 
e.g., Eq. (17) from Mankofsky et al., 1987). 

For massless electrons the equation of motion of the electron 
fluid takes the form of the standard generalized Ohm’s law (e.g. 
Braginskii, 1965) 


E= — Ge X B)- — V Pe -^ 
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where p e = nm e < vf >/3 = n e k B T e , and v' e are the scalar electron 
pressure and the thermal velocity of electrons, and the electron 
current is estimated from Eq. (5). 

The induction equation (Faraday’s law) has a form 


ldB 
c at 


+ V x E = 0. 


( 4 ) 


The total current is given by 
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( 5 ) 


where U s is the bulk velocity of ions of the type s. 

Since we suppose that electron heating due to collisions with 
ions is very small, the electron fluid is considered adiabatic. 
For simplicity, we assume that the total electron pressure may 
be represented as a sum of partial pressures of all electron 
populations 


p e cc 


(Pe n lup + Pe, PI 
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( 6 ) 


where p e and /? eP1 denote electron upwind and pickup betas, 
respectively. Note that p e k = p ek /(B 2 /8n), where k is a population 
of electrons. We also assume here that n e , up = n i>up , n e ,pi = n l PI . 

The neutral atmosphere of Europa serves as a source of new 
ions, mainly by electron impact ionization from corotating 
(or nearly corotating) plasma and also by photoionization. The 
neutral atmospheric molecules also serve as collisional targets for 
charge exchange by corotating ions. The impacting ions consist of 
both upstream torus ions and newly implanted ions which are 
picked up by the motional electric field. 

In the current model, we assume that the background plasma 
contains only the ions with molecular mass/charge of 8 and 16 
corresponding to 0 ++ and S + + , respectively. 

We assume that Europa has a radius R E = 1 560 km. We have 
also adopted a two-species description for the neutral 0 2 


exosphere of exponential form (Shematovich et al., 2005) 

^neutral, k ~ ^atmos, /<exp[ (r r exobase>fe )//i atmos>fc ], (7) 

where n atmosk denotes the maximum value of the neutral density 
extrapolated to the exobase (n atm0Si i = 3 x 10 4 cm -3 ; n atm0Si2 = 
8.5 x 10 7 cm -3 ; r exobase tl ^ 1700 km; r exobase2 ~ 1560 km), and 
index k denotes either non-thermal (/<=1) or thermal (/<= 2) 
species. Here the scale heights h at m0 s,i = 200 km and h at m0 s ,2 = 
30 km. 

The production rate of new ions from the exosphere near 
Europa corresponds to 

Gexo,k °C ^i,k^atmos,/<:^XP[ — (V— r exobase, kV^atmos./d* (3) 

where n atmos /< denotes the value of the neutral component density 
at r = r eX obase,/< and y ik is the effective ionization rate per atom 
or molecule of species k. v kk includes the photoionization v ph , and 
the electron impact ionization by the magnetospheric electrons 
v e ,im- We assume that our model of the atmosphere mainly 
consists of 0 2 , and we use the effective photoionization rate 1.7 x 
10 -8 s -1 (Johnson et al., 2009). We also adopt the effective 
electron impact ionization rates of 2.4 x 10 -8 cm 3 /s (for 20 eV 
electrons) and 1.1 x 10 -7 cm 3 /s (for 250 eV electrons) (see e.g. 
Johnson et al., 2009). Since the hot electrons represent only 5% of 
the total electron density (see Voyager 1 plasma science (PLS) 
measurements analyzed by Sittler and Strobel, 1987; Bagenal, 
1994) we use the same composition for computing the impact 
ionization rate. We assume that the Sun is located in the direction 
opposite the x-axis. 

The interaction of ions with neutral particles by charge 
exchange (see Eqs. (12)-(15) from Lipatov and Combi, 2006) 
currently includes for the following reactions; 

0 ++ +0 2 ^0 + +0 2 + 

S ++ +0 2 ^S + +0+ 

+0 2 — »0 2 + 02^ (9) 

The effective cross section for charge exchange (a c , ex = 2.6 x 
10 -19 m 2 ) was the same as that used in the hybrid modeling of 
Io’s plasma environment (see Lipatov and Combi, 2006; McGrath 
and Johnson, 1989). A more complete list of reactions will be 
considered in future modeling. Of course, this also requires the 
addition of Monte Carlo computations. However, this approach is 
beyond the scope of this paper. 

Our code solves Eqs. (l)-(9). 

We discuss two models of the interaction between the Jovian 
magnetosphere and Europa. In Section 3.1 we discuss the inter- 
action model for the cases with different ions and electron betas, 
different pickup ion production rates near the surface of Europa, 
and homogeneous global Jovian magnetic field (model I, cases 
(a) and (b)), whereas in Section 3.2 we consider the model II, case 
(c) with realistic global Jovian magnetic field and the internal 
dipole magnetic field placed in the center of Europa. To study the 
interaction of the plasma torus with the ionosphere of Europa, the 
following set of Jovian plasma torus and ionosphere parameters 
were adopted in accordance with the Galileo Europa E4 flyby 
observational data (Paterson et al., 1999; Khurana et al., 1998; 
Kivelson et al., 1997, 1999): magnetic field, B 0 = 469nT and 
B = (77.6,-140.7,-441.3) nT; torus plasma speed relative to 
Europa (Paterson et al., 1999), lf 0 = 105km/s; upstream ion 
densities, p 0 ++ =10 cm -3 ; p s ++ =10 cm -3 and ion temperature, 
T( = 25-100 eV (Paterson et al., 1999); electron temperature for 
suprathermal population, T e = 20 eV (Sittler and Strobel, 1987); 
ratio of specific heats, y = 5/3; Alfven and sonic Mach numbers, 
M a = 0.25, M s = 3.66. 
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Initial conditions: Initially, the computational domain contains 
only supersonic and sub-Alfvenic plasma torus flow with a 
homogeneous spatial distribution and a Maxwellian velocity 
distribution; the pickup ions have a weak density and spherical 
spatial distribution. The magnetic and electric fields are B = B 0 
and E = -U 0 x B 0 . Inside Europa the electromagnetic fields are 
E = 0 and B = B 0 , and the bulk velocities of ions and electrons are 
also equal to zero. Here the x-axis is directed in the corotation 
direction, the y-axis is directed toward Jupiter, and the z-axis is 
directed to the north, as shown in Fig. 1. In model I, cases (a) and 

(b) we use a homogeneous magnetic field for the initial and 
boundary conditions (see paragraph above). In model II, case 

(c) we use an extrapolation of the magnetic field profile along E4 
trajectory (see, Kivelson et al., 1999, 2009) onto the computation 
domain for the initial and boundary conditions. The effect of 
global variation in the magnetic field in the rest of Europa was not 
taken into account directly in the modeling but it was included in 
to the modeling as an internal magnetic dipole (see, Schilling 
et al., 2007, 2008). 

At t > 0 we begin to inject the pickup ions with a spatial 
distribution according to Eq. (8). Far upstream (x = -15 R E ), the 
background ion flux is assumed to have a Maxwellian velocity 
distribution. 

Boundary conditions: On the side boundaries (y= ±DY/2 and 
z= + DZ/2), periodic boundary conditions were imposed for 
incoming flow particles. The pickup ions exit the computa- 
tional domain when they intersect the side boundary surfaces 
y = DY/ 2-5 x Ay, y=-DY/ 2 + 5 x Ay, z = DZ/2- x A z, z=-DZ/ 
2 + 5 x Az. Thus there is no influx of pickup ions at the side 
boundaries. 

At the side boundaries, we also use a damping boundary 
condition for the electromagnetic field (see e.g., Lipatov and 
Combi, 2006; Umeda et al., 2001). This procedure allows us to 
reduce outcoming electromagnetic perturbations, which may be 
reflected at the boundaries. 

Far downstream (; x = 12 R E ), we adopted a free escape condition 
for particles and the “Sommerfeld” radiation condition for the 
magnetic field (see e.g., Tikhonov and Samarskii, 1963) and a free 
escape condition for particles with re-entry of a portion of the 
particles from the outflow plasma. 

At Europa’s surface, r = R E « 1560 km, the particles are 
absorbed. In model I, there is no boundary condition at Europa’s 
surface for the electromagnetic field; we also use our equations 
for the electromagnetic field, (see, Eqs. (2), (4) and (9) from 
Lipatov and Combi, 2006) inside Europa but using the low 
internal conductivity (Reynolds number, Re = 0.5) and very small 
value for bulk velocity that are calculated from the particles. 
In model II, we also use an inductive magnetic dipole 


(0,0, -72.5) nTR^ for the boundary condition at Europa’s surface 
that simulates the effect of nonstationarity of Jovian magnetic 
field at the position of Europa. In this way, the jump in the electric 
field is due to the variation of the value of the conductivity and 
bulk velocity across Europa’s surface (note that the center of 
Europa is at x = 0 ,y = 0,z = 0). 

The three-dimensional computational domain has dimensions 
DX = 27R e , DY = 30 R e and DZ = 30R E . We used mesh of 
301 x 301 x 271 grid points, and 5 x 10 8 and 5 x 10 8 particles 
for ions and pickup ions, respectively, for a homogeneous mesh 
computation. The particle time step A t p and the electromagnetic 
field time step A t EB satisfy the following condition; v max A t v < 
min(Ax,Ay,Az)/8 and v max At E B < min(Ax,Ay,Az)/256. 

The global physics in Europa’s environment is controlled by a 
set of dimensionless independent parameters such as 1W A , /? e , 
Mi/Mp, ion production and charge exchange rates, diffusion lengths, 
and the ion gyroradius e = p d /R E . Here p ci = H 0 / ( eB/M[C) = 
M A c/co pi and the ion plasma frequency co P i = ^/47inoe 2 /M[. Mi 
and M p denote the ion and proton masses. For real values of the 
magnetic field, the value of the ion gyroradius is about 80 km, 
which is calculated from the local bulk velocity. The dimension- 
less ion gyroradius and grid spacing have the values e = 0.05 and 

a x /r e = 0.1. 

In order to study ion kinetic effects (e.g. excitation of low- 
frequency oscillations (co < Q b ) by the mass loading), we must 
satisfy the condition A < (10-20 )c/co pb , where Q b and co pb denote 
the gyrofrequency and the plasma frequency for upstream ions 
(Winske et al., 1985). The above estimation of the plasma 
parameters shows that we have good resolution for the low- 
frequency waves (see also Lipatov et al., 2012). 

There is another problem — numerical resolution of the gyro- 
radius on the spatial grid. It becomes very important near 
Europa’s surface where the MHD model cannot to be used and 
we have to use a kinetic model to study a trajectory of the heavy 
ions and their interaction with the surface of Europa. Our current 
model still does resolved the last effect and we expect to 
improved the model by use a spherical system of coordinates in 
future research. 


3. Results of Europa’s environment simulation 

3.1. Effects of plasma betas on the plasma wake structure 

In order to study the effect of plasma parameters on the 
structure of the plasma wake and the Alfven wing, we have 
performed a modeling (model I) for two cases (a) and (b) with 
different values of the upstream background ion temperatures, 
pickup electron temperatures, and value of the pickup production 
rate near the surface of Europa. 

The following plasma parameters are chosen the same for both 
models; full magnetosphere corotation speed is U 0 = 105km/s; 
upstream densities are p 0 ++ =10 cm -3 , p s ++ =10 cm -3 ; mag- 
netic field is B 0 = 469 nT; B = 77.6,-140.7,-441.3 nT; Alfvenic 
Mach number M A =0.25; magnetosonic Mach number M s = 3.66. 
The model of 0 2 atmosphere was taken from Cassidy et al. (2007), 
Shematovich et al. (2005) and Smyth and Marconi (2006). In 
model I, cases (a) and (b), Europa’s interior is represented as low 
conducting body with Reynolds number Re = 0.5. 

Model I , case (a): Upstream ion temperatures are T 0 ++ = 
25 eV; T s ++=25eV and upstream electron temperature is 
T e> o = 20 eV. Temperatures of electrons connected with non- 
thermal and thermal pickup ions are T ejlon . tberma i = 20 eV; 
Te, thermal — 20 eV. 
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X/Re Z/Re X/Re 


Fig. 2. 2D cuts of non-thermal pickup ion 0^ density profile. Model I, case (a) (top) and case (b) (bottom), x-y cuts (left column) are located at z= 0, y-z cuts are located 
at x/Re = 7, and x-z cuts (right column) are located at y=0. 


a 


N(02+, thermal) N(02+, thermal) 


N(02+, thermal) 



-15 -10 -5 0 5 10 -10 -5 0 5 10 15 -15 -10 -5 0 5 10 

X/Re Z/Re X/Re 


Fig. 3. 2D cuts of the thermal pickup ion density profile. Model I, case (a) (top) and case (b) (bottom), x-y cuts (left column) are located at z=0, y-z cuts are located 
at x/Re = 7, and x-z cuts (right column) are located at y=0. 


Model I , case (b) (reduced density for thermal 0 2 by factor 60 
near surface and higher electron temperatures; increased upstream 
ion temperatures, T o ++=100eV; T s ++=100eV): upstream 


electron temperature is T e>0 = 20 eV; temperatures of electrons 
connected with non-thermal and thermal Oj pickup ions 

T e, non-thermal — 200 eV, T e, thermal — 200 eV. 
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a N(0++) N(0++) N(0++) 



X/L Z/L X/L 


Fig. 4. 2D cuts of the background 0 + + ion density profiles. Model I, case (a) (top) and case (b) (bottom), x-y cuts (left column) are located at z=0,y-z cuts are located at 
x/R e = 7, and x-z cuts (right column) are located at y=0. 


We have performed a several hybrid modeling with different 
ions and electron betas, and different rate productions for 
pickup ions but we discuss here only the modeling, which can fit 
the observations. 

The initial thermal velocities of 0^ non-thermal and thermal 
ions are chosen as the following: v thtTlon _ therma i = 3.0 km/s (2 eV) 
and v th dermal = 0.5 km/s (0.05 eV). The initial bulk velocity of 
pickup ions is about 1 km/s. Eq. (8) gives the following total 
pickup ion production rate: Qo+, thermal — 0.825 x 10 28 s -1 and 
Q-O* , non-thermal = 1 .95 X 1 0 26 S” 1 . 

Let us consider first the global picture of the interaction of the 
plasma torus with Europa. The results of this modeling are shown 
in Figs. 2-4. Figs. 2 and 3 demonstrate 2D cuts for non-thermal 
and thermal pickup ion density profiles. One can observe the 
asymmetrical distribution of the pickup ion density (top, case (a)) 
and (bottom, case (b)) in the x-y, y-z (x = 5R E ) and z-x planes. 
The pickup ion motion is determined mainly by the electromag- 
netic drift. The motion along the magnetic field is due to the 
thermal velocity and the gradient of the electron pressure. A more 
wider density profile of the pickup ions was observed in the case 
(b), Figs. 2 and 3 (bottom). 

The figures demonstrate a strong structuring in the non- 
thermal and thermal 0^ ion density profiles. While case 
(a) produces a much higher peak in the thermal 0^ ion density 
as was seen in E4 observations, case (b) produces much better 
agreement with observation for the thermal 0 } ion density as 
shown in Figs. 2 and 3. 

The modeling also demonstrates the asymmetrical distribution 
of the background 0 + + ion density in the x-y , y-z (x = 5 R E ) and 
z-x planes, Fig. 4. The asymmetrical distribution of the back- 
ground ions in the x-y plane may be explained by the existence of 
a strong B z component in the upstream magnetic field. One can 
also see an increase in the plasma density near Europa due to the 
formation of a magnetic barrier (not shown here). In case (b) this 




y=0, z= 0. Model I, case (a) (top) and case (b) (bottom). 

effect is stronger than in case (a). The density profiles for SO + + 
background ions are close to the density profiles for 0 ++ ions. 

The inclination of the magnetic field results in an asymme- 
trical boundary condition for ion dynamics (penetration and 
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Y(Re) 





Fig. 6. ID cuts of the background 0 + + /S + + , and pickup non-thermal/thermal (0^) ion densities from simulation. Y(Re ) denotes a projection of satellite trajectory onto 
the y axis. Model I, case (a) (left) and case (b) (right). Bottom — E4 observation of the total ion density (Paterson et al., 1999). 


reflection) in Europa’s ionosphere and an asymmetrical Alfven 
wing. 

Note that the background ion flow around the effective 
obstacle that is produced by pickup ions and the ionosphere. 
The pickup ions flow from the “corona” across the magnetic field 
due to electromagnetic drift, whereas the motion along the 
magnetic field is determined by the thermal velocity of ions and 
the electron pressure. 

Fig. 5 demonstrates the ID cuts (y=0, z= 0) of the background 
density 0 ++ for case (a) (top) and case (b) (bottom). Strong 
jumps in the plasma density with N 0 ++ max = 80 cm -3 (case (a)) 
and N 0++max = 17 cm -3 (case (b)) are observed on the day-side of 
the ionosphere, whereas a reduction in the plasma density is 
observed in the plasma wake. Note that the jump in the plasma 
density profile is stronger in case (a) than it is observed in case (b). 
Both jumps are located near the surface of Europa. 

Fig. 6 shows ID density profiles of the background and pickup 
ions along the E4 trajectory of the Galileo spacecraft. One can see 
a strong plasma void in the center of the plasma wake. There is 
also a sharp boundary with an overshoot in the density profiles on 
the side of the plasma wake in the Jupiter-direction, and a smooth 
boundary layer on the side in the anti-Jupiter direction, Fig. 6 
(top). The density profile for 0 + + is similar the density profiles 
for the S + + upstream ions. Fig. 6 (middle and bottom) also shows 
the density profiles for the non-thermal (top) and thermal 
(bottom) Oj pickup ions. One can see the split structure of the 


plasma tail. The effect of splitting of the plasma tail was also 
observed in the hybrid modeling of weak comets (see, e.g., Lipatov 
et al., 1997; Lipatov, 2002). The general feature of this plasma 
density is due to the effect of the finite heavy gyroradius. The total 
ion density profile observed in E4 pass is shown in Fig. 6 (bottom). 
The observed value of the density in these peaks is lower than in 
modeling and it may be explained by an overestimated density of 
0^ pickup ions for case (a). In the case (b), disagreement is not as 
strong, an improvement of the atmosphere model is still required. 

The modeling gives the following total fluxes for the 
pickup ions (case (a)): 1.4 x 10 22 mol/s (non-thermal) and 1.75 x 
10 25 mol/s (thermal); (case (b)); 0.8 x 10 22 mol/s (non-thermal) 
and 1.0 x 10 25 mol/s (thermal) across the back boundary 
x=UR e . 

Let us consider a global distribution of the electric and 
magnetic field in Europa’s environment. Fig. 7 shows B x , B z 
magnetic and E y electric field profiles for case (a) (left) and case 
(b) (right). The y-z cuts (top and middle) are located at x/R E = 7, 
and x-y cuts (bottom) are located at y=0. The figure demon- 
strates perturbations in the magnetic B x and electric E y field 
profiles, which are due to the formation of an Alfven wing. 
The increase in the magnetic field B z indicates the formation of 
an asymmetrical magnetic barrier, Fig. 7 (bottom). 

The asymmetry of the modeling distributions in B appears to 
be caused by the finite gyroradius effects of incoming and pickup 
ions. A weak perturbation of the magnetic field was observed near 
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Fig. 7. 2D cuts of the B x , B z magnetic and E y electric field profiles. Model I, case (a) (left) and case (b) (right), y-z cuts (top and middle) are located at x/R E = 7, and y-x cuts 
(bottom) are located at z/R E = 0. 


the ionosphere of Europa: compression of the upstream magnetic 
field and decompression in the plasma wake. 

The modeling also shows the formation of an Alfven wing in 
the direction of the main magnetic field. The formation of the 
Alfven wing in a sub-Alfvenic flow near Europa is similar to a 
formation near Io, which was first studied analytically by 
Neubauer (1980). The pickup ions play important role on the fine 
structure of the Alfven wing due to effects of mass loading. In 
particular, the scale of the front of the Alfven wing must be 
determined by the gyroradius of pickup ions. Unfortunately, in 
our 3D hybrid kinetic simulation we cannot yet resolved this 
spatial scales. 


3.2. Effects of inductive Europa’s magnetic field 

In the first set of models (Section 3.1, model I, cases (a) and 
(b)), we used a homogeneous global magnetic field as an initial 
condition. These models do not produce agreement between the 
simulated and observed magnetic fields. 

In the second set of modeling, we take into account the 
gradient of the global Jovian magnetic field for an initial magnetic 
field distribution. In the paper by Kivelson et al. (1999, 1997, 
2000), it has been shown that the B y component of the magneto- 
spheric magnetic field has strong time variations at the position 
of Europa. In the MHD-fluid approximation, the effects of such 
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magnetic field variations are estimated in Schilling et al. (2007, 
2008). The initial plasma density and bulk velocity distribution 
in our modeling were taken from the E4 flyby data (Paterson 
etal., 1999). 

We created the following model II, case (c), for simulation: the 
density for thermal 0 2 is the same as for model I, case (b), and the 
pickup electron temperature is lower than in model I, case (b). 
The plasma density and bulk velocity distribution in our modeling 
were taken from the E4 flyby data (Paterson et al., 1999): 
full magnetosphere corotation speed U 0 = 105km/s; upstream 
densities are p 0 + + =10 cm -3 ; p s + + =10 cm -3 ; upstream ion and 
electron temperatures, T o ++=100eV; T s ++=100eV; T e >0 = 
20 eV. The temperatures of electrons connected with non- 
thermal and thermal pickup ions are T e non . t h errna i = 100 eV; 
Te, thermal = 100 eV. 

In our hybrid kinetic modeling (model II) we use a simple 
magnetic dipole model of the induced oceanic magnetic field 
from the 10-hour corotation variation of the background Jovian 
magnetic field at Europa (see paragraph “Boundary Conditions”, 


Section 2). And, finally, we fit the results of modeling to the 
components of the measured magnetic field. 

This is not yet a fully self-consistent approach but provides a 
first approximation. Also, the ocean may not be exactly a spherically 
symmetric conducting shell and may ultimately require a higher- 
order multipole model for the induced fields. 

Fig. 8 demonstrates the 2D cuts for non-thermal and thermal 
02 “ pickup ion densities. The figure does not show any extension 
of the pickup ion profile in the y and z directions. The plasma 
wake is narrower in y and z directions in compare with that was 
produced by model I, cases (a) and (b). The reason of this effect is 
due to lower temperature of electrons connected with pickup O^ 
ions than it was in a case (b) and lower pickup ions production 
rate near the surface of Europa than it was in a case (a). 

Fig. 9 shows the distribution of the 0 ++ ion density in the x-y , 
y-z (x = 5 Re) and z-x planes. The narrow plasma wake may be 
explained by smaller temperature of the electrons connected with 
0^ pickup ions, and, hence, with a smaller polarization electric field 
which is responsible for an expansion of Europa’s ionosphere. 
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Fig. 8. 2D cuts of the pickup ion density profile. Non-thermal (top), thermal (bottom). Model II, case (c). x-y cuts (left column) are located at z= 0, y-z cuts are 
located at x/Re = 7, and x-z cuts (right column) are located at y=0. 
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Fig. 9. 2D cuts of the background 0 + + ion density profiles. Model II, case (c). x-y cuts (left column) are located at z=0, y-z cuts are located at x/R E = 7, and x-z cuts 
(right column) are located at y=0. 
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One can also see an increase in the plasma density near Europa 
due to the formation of a magnetic barrier (not shown here). 
The density profile for SO ++ background ions is close to the 
density profile for 0 ++ ions as in model I, cases (a) and (b). 

Fig. 10 shows a ID cut of the background 0 ++ density along 
the x-axis (y = 0,z = 0). One can see jump in the background 
plasma density with N 0 ++ max = 90 cm -3 (model II, case (c)) on 
the day-side of the ionosphere and depletion in the plasma 
density in Europa’s plasma wake. Note that the jump in the 
plasma density profile is stronger in model II, case (c), than it is 
observed in model I, case (a). The jump is located near the surface 
of Europa as it was observed in model I, cases (a) and (b). 

Fig. 11 shows ID density profiles of the background and 
pickup ions along the E4 trajectory of the Galileo spacecraft. 
One can see a strong plasma void in the center of the plasma 
wake. There is also a sharp boundary with an overshoot in the 
density profiles on the left side of the plasma wake, and a smooth 
boundary layer on the right side, Fig. 1 1 (top). The density profile 
for S + + is similar to the density profile for 0 + + background ions. 
Fig. 10 (middle) shows the density profiles for non-thermal and 
thermal 0 2 + pickup ions. The total ion density profile observed 
during the E4 pass is shown in Fig. 1 1 (bottom). Again, one can see 
two peaks in the total ion density profile. However, the observed 
value of the density in these peaks is lower than predicted by the 
model; this may be explained by an overestimated density of 0^ 
pickup ions for model II, case (c). 

The modeling shows that the shaping of Europa’s global 
plasma environment depends on a combination of the upstream 
plasma parameters and pickup ions and electron parameters. For 
example, the reducing in the temperature of electrons connected 
with pickup ions results in the higher density of the thermal 0^ 
pickup ions at the trajectory of a spacecraft, compare Fig. 6 (right) 
and Fig. 11. This effect is connected with the polarization electric 
field which is proportional to the gradient of the electron 
pressure. The reducing in the temperature of the background 
upstream ions results in the widening of the plasma wake, 
compare Fig. 6 (left and right, top) and Fig. 11 (top). These effects 
were earlier demonstrated in the 3D hybrid simulation of Io’s 
plasma environment (Lipatov and Combi, 2006). We have found 
the similarities between the plasma environments of these 
objects. Indeed, Io and Europa have sufficiently thin exospheres 
and strong magnetic fields resulting in a small value of the ion 
gyroradius. 

Let us consider the global distribution for the electromagnetic 
field of model II, case (c). Fig. 12 shows 2D cuts for the magnetic 
B x , B z and electric E y field profiles. The distributions for the B Zl 
E y field shown in the figure are close to the distributions for 
model I, case (b). However, there are significant differences 
between the B x profiles for model I, case (a) and case (b), and 
model II, case (c). The differences between the B x profiles for cases 



y= 0, z= 0. Model II, case (c). 


0++(solid line), S++(bullet) 



total ion density (E4) 



Fig. 11. Background (S + + , 0 + + ), non-thermal (0^ ) and thermal (0^ ) pickup ion 
densities from simulation. Y(l?e) denotes a projection of the spacecraft position 
onto the y axis. Model II, case (c). Bottom — E4 observation of the total ion density 
(Paterson et al., 1999). 

(a) and (b), Fig. 7 (top) are due to much higher density of the 
thermal 0 2 + pickup ions in the plasma wake, whereas the 
differences between the B x profiles for cases (b) and (c) are due 
to the nonlinear interaction of the Alfven wing with nonhomo- 
geneous Jovian magnetic field in model II, case (c). 

Fig. 13 shows the magnetic field components (solid line) B x , B y , 
B z , and \B\ along the E4 trajectory of the Galileo spacecraft. The 
magnetic field components of the inductive magnetic dipole that 
simulates the effect of the nonstationarity of the Jovian magnetic 

field are shown by a dotted line ( ). The circles (o) denote 

observational data from Kivelson et al. (1997) and the initial 
Jovian magnetospheric field at the position of Europa (+++ ). The 
simulation produces a satisfactory agreement with the observa- 
tional data for the B y magnetic field component, but not for the 
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Fig. 12. Model II, case (c). 2D cuts of the magnetic B x (top), electric E y (top) and magnetic B z (bottom) field profiles, y-z cuts (top) are located at x/R E = 7, and y-x cuts 
(bottom) are located at z/R E = 0. 


B x and B z magnetic field components. A multipole model for the 
oceanic magnetic field may address this issue. We will need to 
improve the model of the 0 2 atmosphere, the resolution of the ion 
trajectory, and the gradient in the atmosphere/ionosphere density 
profiles near the surface of Europa to obtain better agreement in 
the B x and B z magnetic field components 


4. Conclusions 

Hybrid modeling of Europa’s plasma environment for E4 
encounter with three ion species demonstrated several features: 


• The modeling shows a strong phase mixing in the plasma 
wake. The plasma wake demonstrates the formation of time- 
dependent structuring in the pickup ion tails (see, e.g., 
McKenzie et al., 2001 for weak comet case) and the splitting 
of the pickup ion tails. The splitting of the plasma wake has the 
same nature as the splitting of the weak comet’s plasma wake 
or the splitting of Titan’s plasma wake. Such finite gyroradius 
effects were also observed in 2.5D hybrid and bi-fluid 


modeling of a weak comet (see, e.g., Lipatov et al., 1997; Sauer 
et al., 1996, 1997; Lipatov, 2002) and in 3D hybrid modeling of 
Titan’s plasma environment (Lipatov et al., 2011, 2012). The 
further investigation of these fine structure needs an addi- 
tional modeling with much better resolution. 

• The model shows the magnetic field barrier formation at the 
day-side portion of the ionosphere. The formation of an Alfven 
wing in the plane of the external magnetic field was also 
observed. Note that the Alfven wing was earlier observed in a 
hybrid simulation of the plasma environment of Io and Europa 
by Lipatov and Combi (2006) and by Lipatov et al. (2010). An 
MHD simulation of the plasma environment of Io and Europa 
also produces the formation of an Alfven wing (Saur et al., 
1999, 1998; Liu et al., 2000; Schilling et al., 2008). 

• The ion and electron temperatures play an important role in 
plasma structure formation, and in creating the ion fluxes 
inside the ionosphere. These effects were observed earlier in a 
3D hybrid simulation of Io’s plasma environment (Lipatov and 
Combi, 2006). Hybrid model produces a correct pickup ion flux 
along the magnetic field in opposite the MHD models which 
operate with pickup ions with a Maxwellian velocity distribu- 
tion. In the current paper, we have presented only three runs 
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X(R e) 

Fig. 13. Magnetic field component profiles along the E4 trajectory after fitting 
with inductive dipole magnetic field. Solid line — modeling, (-) denotes dipole 
field, and ( + ) is the Jovian magnetic field at the position of Europa. o — Galileo’s 
E4 flyby measurements (Kivelson et al., 1997). X(Re ) denotes a projection of the 
spacecraft position onto the x-axis. 

with different combinations of the upstream ion and pickup 
electron temperatures. 

• The model’s total ion density in the plasma wake does not 
satisfactory match the observed density. 

• The constant induced dipole moment (model II, case (c)) 
improves a fit the magnetic field B y component to the E4 
trajectory. However, a fit for the magnetic field B x component 
is still not satisfactory due to the not perfect model of the 
atmosphere/ionosphere and not satisfactory numerical resol- 
ving of the gyriradius on the grid cell. 

• Use of an inhomogeneous background magnetic field provides 
a good agreement between the observed and simulated 
magnetic fields. However, we still need to improve the 
resolution of the gradient in the atmosphere density, the 
gyroradius of pickup ions, and the fields in the internal non- 
conduction ice shell and conduction ocean layers of Europa. 

In our future computational models, we plan to include a 
nonstationary boundary condition for the magnetic field in order 
to take into account the spatially inhomogeneous and nonsta- 
tionary background Jovian magnetic field. This model will also 


respect to a potentially nonspherical ocean shell. We also plan the 
use of a varying atmospheric density, a varying electron tem- 
perature (that plays key-role in the pickup ion dynamics), and 
sputtering processes (Johnson, 1990; Johnson et al., 1998) at the 
surface of Europa. We also plan to use a composite grid structure 
using the “cubed sphere” technique (see, e.g. Koldoba et al., 2002) 
to improve the resolution of the a small scales near the surface of 
Europa and to increase the size of the computational domain. 

The composite grid structure will allow us to estimate the 
inductive magnetic field from the ocean as a part of the total 
current closure that also includes the external plasma currents. 
This technique will allow us to study wave-particle interaction 
effects in the far plasma wake, such as ion cyclotron waves that 
have been observed in the Galileo flyby mission (see e.g. Volwerk 
et al., 2001 ; Kivelson et al., 2009). These models must include the 
induced magnetic field from a putative subsurface ocean, and will 
also include particle trajectory tracing for test particles, e.g. 
electrons and high-energy ions. 

Note that the larger computational domain allows us to use 
the upstream parameters for the plasma and electromagnetic 
field instead of the use of the “damping” boundary condition. 
However, in the outer region of the computational domain (large 
cell size) we have to use a drift-kinetic approach (see e.g. Lipatov 
et al., 2005) for ion dynamics since we cannot approximate the 
ion trajectory there. We can also use a complex particle kinetic 
technique (see e.g. Lipatov, 2012) which provides a flexible fluid/ 
kinetic description and may significantly save computational 
resources. 
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